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This  award  has  allowed  to  develop  a  new,  mathematically  sound  algorithm  of  the  blind  image 
deconvolution,  which  demonstrates  a  superior  performance  when  compared  with  existing  techniques. 

The  concept  change  in  our  approach  to  blind  image  deconvolution  is  in  the  basic  strategy  of  addressing 
sensible  approximate  solutions  to  the  ill-posed  nonlinear  inverse  problem.  These  solutions  are  addresses 
as  fixed  points  of  the  iteration  which  consists  in  alternating  approximations  (AA)  for  the  object  and  for  the 
PSF  performed  with  a  prescribed  number  of  inner  iterative  descents  from  trivial  (zero)  initial  guess.  This 
approach  shall  be  contrasted  with  traditional  inexact  alternating  minimization  (AM)  approach,  where  a 
stationary  point  of  the  objective  function  is  being  targeted  by  the  global  descent  trajectory,  and  monotonic 
descent  to  this  point  is  terminated  to  achieve  regularization. 

Artificial  inversions  with  noise-free  images  show  that  the  new  approach  allows  the  successful 
deconvolution  of  data  of  higher  complexity,  where  the  AM  approach  suffers  from  stagnation.  In  our  vision, 
the  stagnation  is  caused  by  approaching  a  local  critical  point  of  the  cost  function;  these  points  are  not 
necessarily  the  fixed  points  of  our  iteration.  Inversions  with  real  data  produce  solutions  which  are  more 
accurate,  free  from  the  artifacts,  and  do  not  depend  on  the  initial  guess. 

Below  is  a  brief  overview  of  the  algorithm.  Its  complete  description  can  be  found  in  S.  V.  Vorontsov  and  S. 
M.  Jefferies  2015  "A  new  approach  to  blind  deconvolution  of  astronomical  images".  This  paper  has  just 
been  submitted  for  publication  to  the  Inverse  Problems ,  and  is  attached  to  this  Report. 

THE  ALGORITHM  {Alternating  Approximation  method  for  blind  deconvolution  x*y=f} 

Choose  initial  pair  of  regularization  parameters  (kx,ky) 

Choose  non-zero  initial  guess  for  point-spread  function  y 
REPEAT  {Graduated  optimization} 

REPEAT  {Fixed  point  iteration} 
x=x(f;y)  using  kx  descents  with  +SOR  from  zero  guess 

y=y(f;x)  using  ky  descents  with  +SOR  from  zero  guess 

UNTI  L  convergence 
Increase  (kx,ky) 

UNTI  L  stopping  condition  is  met. 

Inner  minimization  routine  performs  standard  (non-blind)  deconvolution,  using  non-negativity- 
constrained  successive  over-relaxation  technique  (+SOR)  with  adaptive  scheduling  of  the  relaxation 
parameter.  This  routine  was  inherited  from  our  pre-award  research,  and  was  added  with  some  few 
improvements  and  with  deeper  convergence  analysis.  In  our  numerical  experience,  this  routine  is 
exceptionally  efficient  and  robust  when  working  with  data  of  arbitrary  complexity.  A  particularly  attractive 
feature  of  the  routine  is  that  it  is  free  from  any  tunable  parameters — apart,  of  coarse,  from  the  number  of 
iterations  where  the  descent  is  terminated  to  prevent  excessive  noise  magnification.  A  natural  stopping 
condition  is  the  flattening  of  the  Fourier  amplitude  of  the  residuals  to  near  the  noise  level. 

Fixed-point  iteration.  The  result  of  the  non-blind  deconvolution  for  object  x  with  given  PSF  y  using  kx 
iterations  of  +SOR,  x=x(f;y)  is  used  in  non-blind  deconvolution  for  y  with  ky  descents,  y=y(f;x),  and  the 

process  is  repeated  until  convergence.  In  our  algorithm,  x  belongs  to  a  compact  convex  set,  and  the 
existence  of  a  fixed  point  of  function  x(f;y(f;x))  is  warranted  by  the  Brouwer  fixed  point  theorem.  We 
consider  attractive  fixed  points  corresponding  to  different  values  of  kx  and  ky  as  candidates  for  a 

sensible  solution  to  the  blind-deconvolution  problem.  The  solution  becomes  independent  on  the  initial 
guess  when  the  first  attractor  is  reached.  To  make  a  distinction  with  AM  methods,  our  approach  to  blind 
deconvolution  in  its  general  form  (+SOR  can  be  replaced  by  another  non-blind  deconvolution  routine) 
shall  be  probably  designated  as  an  "alternating  approximation"  (AA)  method.  Interestingly,  we  have  not 
met  any  reports  on  implementing  this  approach  in  literature  on  nonnegative  matrix  factorization  -  a  wide 
area  of  research  with  numerous  applications,  our  blind-deconvolution  problem  being  just  a  particular 
example. 
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Graduated  optimization.  For  a  given  pair  of  regularization  parameters  (kx,ky),  the  iteration  x<-x(f;y(f;x)) 

may  have  multiple  attractive  fixed  points.  Our  numerical  experiments  show  that  multiplicity  of  the 
attractors  grows  when  regularization  is  relaxed  (bigger  kx,ky),  and  when  addressing  data  of  higher 

complexity  (e.g.  with  a  PSF  of  more  complicated  morphology).  To  find  those  attractors  which  correspond 
to  sensible  approximate  solutions,  we  start  with  an  easier  problem  of  finding  an  appropriate  attractor  at 
small  values  of  (kx,ky),  and  then  relax  the  regularization  gradually  by  going  to  bigger  values  of  (kx,ky), 

using  the  PSF  of  the  previous  fixed  point  as  the  initial  guess.  This  approach,  which  is  known  in  the 
literature  as  graduated  optimization,  is  closely  related  with  approach  known  as  deterministic  annealing. 

When  working  with  data  of  moderately  small  complexity,  this  approach  demonstrates  excellent  results. 
Like  with  any  other  methods,  it  has  to  reveal  its  own  problems  and  limitations  when  the  data  complexity 
becomes  high.  We  have  addressed  these  limitations  by  performing  numerous  artificial  inversions.  When 
image  complexity  is  characterized  by  the  parameter  D/rg,  where  D  is  the  aperture  size  and  rg  is  a 

coherence  scale  of  atmospheric  turbulence,  the  technique  allows  decent  deconvolution  of  single  images 
with  D/rg  up  to  about  20.  This  limitation  is  apparently  related  with  either  some  imperfections  of  the  inner 

minimization  routine,  or  with  too  wide  steps  in  the  regularization  parameter  space  (integer  kx,ky).  Further 

enhancement  of  the  performance  of  the  technique  needs  more  efforts.  We  note,  however,  that  these 
efforts  are  only  relevant  to  practical  inversions  when  data  has  an  exceptionally  high  signal-to-noise  ratio. 

Numerical  results.  We  were  focused  largely  on  the  most  difficult  problem  of  single-frame  blind 
deconvolution,  which  is  important  for  data  sets  that  have  objects  that  change  pose  quickly  and  for  which 
multi-frame  approaches  are  not  valid. 


This  figure  illustrates  one  of  the  results  of  our  artificial  experiments  with  simulated  data  (here, 
D/rg  =10),  targeted  at  examining  the  potential  capability  of  the  deconvolution  technique  to 
accurately  split  the  information,  contained  in  the  blurred  images,  between  the  object  and  the  PSF. 
More  results  can  be  found  in  our  paper.  There,  we  also  report  the  results  obtained  with  two  sets 
of  satellite  images  acquired  using  ground-based  telescopes  with  and  without  adaptive  optics 
compensation  at  the  Air  Force’s  AMOS  site  on  Maui,  and  extend  the  technique  to  multi-frame 
deconvolution.  The  results  are  compared  with  those  obtained  by  alternating  minimization  (AM) 
using  the  positivity-constrained  conjugate-gradient  version  of  the  PCID  algorithm,  which  is  in 
routine  use  at  the  AMOS  facility  on  Maui,  demonstrating  much  better  performance  of  the  new 
technique. 
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Further  efforts.  There  are  at  least  three  directions  where  the  technique  calls  for  extensions. 
First,  the  algorithm  in  its  current  form  implements  least-squares  (LS)  approximations,  which 
maximize  the  solution  likelihood  only  when  the  noise  statistics  is  stationary  Gaussian.  In  reality, 
noise  in  the  images  is  close  to  Poissonian  (photon  shot  noise).  Despite  the  observations  that 
accounting  for  the  correct  noise  statistics  does  not  change  the  results  of  image  deconvolution 
significantly  when  compared  with  the  LS  approach,  adaptation  of  the  technique  to  proper 
statistics  has  to  be  done,  as  otherwise  the  results  can  not  be  of  minimum  variance  because  they 
are  not  based  on  a  correct  data  weighting.  The  main  difficulty  is  in  modifying  the  theoretical 
recipe  for  scheduling  the  relaxation  parameter  of  the  +SOR  routine,  which  is  used  in 
minimization. 

When  applied  to  multi-frame  deconvolution,  our  technique  (together  with  others)  attributes 
equal  weights  to  different  data  frames,  and  treats  them  on  an  equal  basis.  Adaptation  has  to  be 
made  to  account  for  the  fact  that  separate  frames  are  often  distorted  by  the  atmospheric 
turbulence  to  a  very  different  extent.  This  problem  is  related  with  a  similar  one  in  the  single- 
frame  deconvolution,  where  an  optimal  relation  between  the  number  of  descents  in  the  object- 
and  in  the  PSF  space  has  to  be  addressed  from  the  theoretical  standpoint. 

The  band-limited  nature  of  astronomical  imagery  (which  comes  from  the  finite  aperture  size) 
has  to  exploited  by  the  algorithm.  This  is  particularly  important  for  attempts  of  improving  the 
resolution  of  distant  objects  to  beyond  the  diffraction  limit  by  using  over-sampling  (collecting 
data  on  a  grid  finer  than  the  diffraction  limit).  An  obvious  way  is  re-parameterization  of  the  PSF 
in  terms  of  the  wave- front  distortion  at  the  aperture  plane.  The  difficulty  to  overcome  is  in  the 
additional  source  of  nonlinearity  introduced  to  the  inverse  problem. 


S.V.Vorontsov 
7  June  2015 
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